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Materials and Methods 1: Derivation and solution of the self-consistent 
equations - Details of the algorithm 

Derivation of the self-consistent multimode laser equations. The Green function which 
forms the kernel of the fundamental self-consistent equation Eq. 1 has a non-hermitian spectral 
representation of the form (SI) 

in=l ^ \ // 

where the CF states (pmix, k) satisfy the equation 
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VVm(a;, k) = ki{k) ip^ix, k) (S2) 



with the boundary condition that outside the gain region there are only outgoing waves with 
wavevector k. The biorthogonal partners ipm{x, k) satisfy the complex conjugate equation with 
incoming wave boundary conditions, hence iprn{x, k) = (e{x)ipm{x, k))*. The biorthogonality 
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condition on these functions is J^^.^ dx <^^(aj, k)Lpn{x, k) = 5mn with appropriate normaliza- 
tion. Writing G{x, x'; k) in Eq. 1 in the form Eq. [STl substituting the CF expansion of the 
unknown lasing modes \I/^(cc) = at^ipmix, k^), using the biorthogonality relation and 

truncating to Ncf CF states nearest the atomic line (gain center) yields Eq. 2, with the non- 
linear operator = Tmnik^; a^) where 



7_L - i{k - ka) 2ka{k - km{k)) 

(1 + do{x'))^l{x', k)Vn{x\ k) 

e{x'){l + T..TA^.{x')\^) ■ ^''^ 

We refer to T^nik; a^) also as Tmn{k) or Tmnik; Dq) in the text in order to emphasize the 
dependence that is relevant in the discussion. 

Threshold Matrices: The set of non-linear equations (2) have the following properties. Be- 
low some finite value of the pump Dq only the trivial solution exists: = 0, VyU. As Dq 
is increased, a series of thresholds are reached at which the number of non-trivial solutions in- 
creases by one. If we denote the threshold D^^ for /-mode lasing there exist I solutions {a^, A;^} 
(/i = 1, . . . , /) for a pump parameter Dq such that D^^ < Dq < D^^^\ In each of these in- 
tervals we assume that / solutions to Eq. 2 exist and find them by a method to be described 
below. 

In order to find the first threshold we consider the linear operator Tmn{k] = 0) = T^^^ {k) 
where 

T^ik) = 7," , . , \ / ^ji + dQ{x'))r ix',k)^^{x',k) ^^^^ 



1± - i{k - ka) 2ka{k - k„i{k)) J e{x 

which is obtained by neglecting the J2u '^u\'^u{x')\'^ term in the denominator of Eq. [S3l The 
resulting linear equation associated with (2) has the form 

r(°)(fc)a'^ = (l/DoX (S5) 

which has solutions (as noted in the text) when an eigenvalue of this matrix, Xn\k),n = 
1, . . . , Ncf, is real and has the value 1/Dq (i.e. DqXu^ = 1). Since T^^\k) is independent 
of Dq, so are its eigenvalues. 7'*-°'' (k) is non-hermitian and has complex eigenvalues for general 
values of k. Denoting the largest eigenvalue of this matrix by \ solve Eq. [S5]by tuning 



k until Im 



\'^'>(k = ki) = 0; this determines the threshold d[1^ = l/\f\k = ki). ki is the 

lasing frequency of the first mode at threshold; the corresponding eigenvector of 'T^^\k) gives 
the projection of the lasing mode onto the CF states, a\ at threshold. The "length" of is 
not determined from Eq.[S5]but rises continuously from zero at threshold and is determined by 
the non-linear equation (2) infinitesimally above threshold. The other, smaller eigenvalues of 
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T^^^ (k) define the non-interacting thresholds for other modes with their frequencies determined 
by the same reality condition, but the actual thresholds of all higher modes will differ substan- 
tially from their non-interacting values due to the non-linear term in Eq. 2 which now comes 
into play. The actual lasing frequencies of higher modes have only a weak dependence on Dq 
and differ little from their non-interacting values (see Fig. 2, inset). 

Above the threshold D^j^ we solve the non-linear Eq. 2 by an iterative method to be de- 
scribed below. Assuming we have this solution in hand at each value of Dq we can construct 
the generalized interacting threshold matrix 

p, N _ ^7± 1 [ , / (1 + dQ{x'))ipl,[x' , k)ipn{x', k) 

' ^^-i{k-ka)2ka{k-k^{k))J e{x'){l + T{k,)\m,{xX) 

Here, the Dq dependence of Tml derives from the non-linear dependence of (assumed to 
be determined by the procedure described further below) on Dq. We will alternatively use the 
notation Tmn{k; a^) in what follows. If we vary k at fixed Dq this linear operator will have 
a real eigenvalue X^i\ki) = I/Dq reflecting the existence of the first lasing mode. To find 
the threshold for the second lasing mode we vary k until its second largest eigenvalue satisfies 

Im A2^'*(/c = ^2) = 0. Similarly to the non-interacting case, k2 is the lasing frequency of the 

second mode and its expected threshold is D[f^ = l/\2\k = ^2)- This procedure generalizes 
in the obvious manner to the third and higher thresholds. Our algorithm for calculating the 
multimode lasing states continuously monitors the relevant threshold matrix as Dq is varied to 
determine at each pump power how many lasing modes are turned on. The eigenvalues of these 
threshold matrices at a fixed Dq have a smooth flow in the complex plane and we can always 
associate a particular eigenvalue at an arbitrary k with a particular lasing mode fi (which may 
not be yet turned on). Henceforth fi is assumed to be ordered according to the order of turn-on. 
The eigenvalue flow of the non-interacting threshold matrix is illustrated in Fig. [STl note that 
the initial values of each of the eigenvalues are different, as shown by the star- shaped markers. 

Non-linear Solver: The interacting threshold matrix T^^\k; Dq) provides us with the start- 
ing values of the lasing frequencies k^ and the threshold solution (up to a proportionality 
constant) infinitesimally above threshold for the non-linear lasing equations (2). Equation 2 
appears to be convenient for iterative solution, but it must be further constrained in order for 
this procedure to work. Tmnik; a^) is invariant under global phase rotations e"^a^ and 

(2) only has a unique solution when this overall phase is fixed (the "gauge" is fixed). This phase 
can be fixed in a trial solution, but it will be changed by each iteration of the equation and so the 
correct procedure is to allow the trial lasing frequency to flow under iteration so as to maintain 
the desired global phase of the trial solution. In this manner the interacting lasing frequencies 
can be found above threshold (note that these frequencies are different from the non-interacting 

or threshold values). In practice, we choose the gauge by setting Im a^j^ = 0, where is 
the largest CF component of the eigenvector of the non-interacting threshold matrix T^"-* (k). 
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With this important modification the solutions to Eq. 2 are found by increasing Dq in small 
steps above the first threshold and solving for the fixed point(s) of the equation, the vectors a'^, 
by iteration. For this we use a multi-dimensional root finder based on the Powell hybrid method. 
Convergence depends on the quality of the initial approximation which in turn depends on how 
fine the pump range is discretized. We find that close to the thresholds the rate of convergence 
is in general slower as would be expected for non-linear systems close to a bifurcation. At each 
value of Dq the associated interacting threshold matrix is constructed from the non-zero a'^ 
which have been found and monitored to check if the next lasing mode has reached threshold 
and should be included in the non-linear system. Note that the number of lasing modes is not a 
monotonically increasing function of Dq; we find that lasing modes can "turn-off due to strong 
modal interactions, as described below (black mode in Fig. 2 and Fig.[S2l). 

The eigenvalues of the interacting threshold matrices as a function of Dq are very interesting 
because they show the strong effects of mode competition. If we plot DqX^/;'^ vs. Dq these are 
just straight lines intersecting unity at the non-interacting thresholds; the interacting eigenvalues 
will be sub-linear, leading to much higher thresholds, and some will even be decreasing with 
increasing Dq, indicating modes which are completely suppressed by mode competition and 
might never turn on. This behavior is shown in Fig. [S2] below along with the behavior of the 
lasing frequencies vs. Dq. Strong mode-mode interactions mediated by gain- saturation can 
be studied in detail in Fig. [S2l For instance we observe that the turn-on of the black mode is 
delayed from Dq / Dqc ~ 72 to Dq / Dqc ~ 78 due to interactions mainly with the orange mode, 
which turns on earlier. From Fig. [S2lB) we see that the frequencies of these two modes shift 
closer to each other, which increases the interaction and results finally in the turn-off of the 
black mode at about Dq/Dqc ~ 117, as seen in Fig. 2 and Fig. [S2l A). At this point we observe 
a kink in the intensity of the orange mode. A similar interaction takes place between the green 
and purple modes as described in the main text. 



Materials and Methods 2: Collective contribution to the laser frequency 

Here we provide details leading to Eq. 3 of the main text; in this section, we will measure 
all frequencies from the atomic transition frequency ka to simplify the equations (e.g. — 



kn). The gauge-fixing condition Im 



leads to the following equation for the 



corresponding lasing frequency (we will set - 

^MRe[A>;]-n>;im[A't]) 



1) 



k,. 



(7x + n't)Re[A>t] - {k, - q'i)Im[A^,] (7^ + - {k, - g^)^/. ' 
Here k'^ = Qi — it^i ( > 0) is the CF frequency of the largest contributing CF component, 

- dQ{x'))^>t\x\kM{x\k) 



(S7) 



NcF 
71=1 



dx' - 



e(^')(l + E.r.|v]/,(a;') 



(S8) 
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and = Im[A'l]/ RelA'^]. Equation|S7]is exact, but it is useful to make minor approximations 
in order to get a more easily interpreted result. In our parameter range is typically much less 
than one and so we can replace in the denominator of Eq.[S7]with its first order approximation 
(i.e. the result for cr^ = 0), which is kj^^ = i+^m^^^ • This leads to the result of Eq. 3 

[l-tani^^']^], (S9) 

where tan[0^] = Im[A'^]/ Re[Ai] and we identify the collective contribution to the lasing fre- 
quencies by kl^'^ = — taii[(j)'^]Ki / Qi . The first term, kj^ \ is well-known from single mode 
lasing, it represents the pulling of the cavity frequency Qi towards the atomic line (ka = in 
our current convention). In the fractional finesse limit of a DRL > 7^, it alone would give 
k^ ~ Q'l (7±//«i) < Qi, i.e. a very large pulling of the lasing frequencies towards the atomic 
line center. This effect is seen in Fig. 1. However the second term in Eq. [S9]is the collective 
effect due to all the other CF states. This contribution has no analog in conventional lasers and 
is random in sign, as can be seen in Fig. 1, where some frequencies are pushed towards and 
others away from the atomic line center due to this term. The size of this effect depends on the 
magnitude of = Im[A'l^]/ Re[Ai]. As noted, in our parameter range this quantity is small; 
this is due to a remnant of the biorthogonality relation. Analysis of the quantity A'^ related to 
suggests that at larger values of kaR and with more lasing modes (higher above threshold) 
the quantity can be large and of arbitrary sign. When this is true the approximation leading 
to Eq.[S9]will not be valid, but from Eq.[S7]we see that the lasing frequencies will be dominated 
by the collective effects of all the CF poles and will not be associated with any single CF state 
or passive cavity resonance. The results discussed in Materials and Methods 4 below support 
this conjecture. 

Materials and Methods 3: Numerical calculation of CF modes of the DRL 
and parameters 

In our model for the DRL the region of uniform gain is assumed to be a disk of radius R, on 
which the differential equation [S2l is discretized with a polar mesh {p{i),(j){j)), i = 1, . . . , Np, 
j = 1, . . . , chosen to be finer than the wavelength of the light \ = 27c/k{kis eventually to be 
set to the respective lasing frequencies k^). The resulting eigenvalue equation is Cij ipmih j) = 
k'Lik)(pmii,j) where Cij = -7(^)^l(i),^ij)^ with V^^.^ ^^^^.^ the discretized Laplacian in polar 
coordinates, see Ref. (52). The static dielectric disorder enters the discretized operator at each 
grid point explicitly via e(pj, and can thus be chosen at will; in what follows the dielectric 
function at each grid point takes the values e = ei = (1.2)^ or e = eg = 1 randomly with 
roughly 20% coverage of "nanoparticles". The outgoing CF boundary condition is imposed by 
continuously connecting the solution of Eq.[S2]and its derivative to a superposition of outgoing 
Hankel functions, Hm\kR)e^^^'^, on the boundary of the gain-disk. The resulting boundary 
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conditions are /c-dependent and can be written as 

^{N, + 1,0,) = J2<^MN,)] ( 1 + kS ^S)fu^l ] ^^^^^ 

where 5 = pi+i — and am[</5(Ai'p)] are the discrete angular Fourier coefficients of the solution 
on the last ring v3(A^p, 0j) 

a^[^(iV,)] = — ^^(iV„0,)e-*-^' (Sll) 

i 

This leads to a finite non-hermitian eigenvalue problem depending parametrically on k which 
is solved by customized linear algebra packages for sparse matrices. Eqs. ISIOI and ISl II are 
incorporated by adding a A^^ x A^",^ block to the discretized Laplacian £y which renders it k- 
dependent. The resulting solutions {k^ik)-, ^mix, k)} are used during the iteration procedure 
in the construction of the non-linear operators Tmn{k] a^). 

The rest of the parameters used in the calculations presented are as following: 7^/2 = 1, 
kaR = 30, NcF = 16. 



Materials and Methods 4: Spatial structure of lasing modes 

The amplitudes a'^ determine the spatial structure of the lasing modes via 

NcF 
m=l 

Since these amplitudes evolve continuously from the relevant eigenvectors of the threshold ma- 
trices (see discussion above), it is useful to consider their distribution from analysis of these ma- 
trices. Construction of the non-interacting threshold matrix T^^^ (for uniform pumping) finds 
this matrix to be close to diagonal and its eigenvectors to be localized on single CF states, imply- 
ing that at the first threshold the lasing state is primarily made up of one CF state. The integral in 
Eq.|S4]is not diagonal by biorthogonality due to the factor e{x') in the denominator. However, it 
can be divided into an integral over the region with background dielectric function eg and scat- 
tering centers with e = ei. Adding and subtracting appropriate quantities and using biorthogo- 

nality one finds that T^, cx / dx' ^^-^^^^^f^ = i [5„„ - ^ /^^ dx'^*^{x', k)y,4x', k) , 
where gi denotes the area of the gain region with e = ei. The second term here is small due 
primarily to the fluctuating phases of the CF states, but also due to the smallness of the prefactor 
and the fact that we have taken gi ^ 20% of the entire gain region. When T'-^-' is approximately 
diagonal, the quantity discussed above is small. As a result, at low pump powers for the first 
lasing mode, /c^ ^ kf!^ with a small collective contribution. 



6 



However, if one analyzes the interacting threshold matrix T^^-* {k; Dq) (and higher ones) one 
immediately sees that the presence of lasing modes modifies the denominator in a crucial man- 
ner: there is now a space-dependent term (the "hole-burning term") which cannot be divided up 
into two constant regions, but instead varies continuously throughout the gain medium. As this 
term increases, the threshold matrices and the non-linear operator T^^^ become more and more 
non-diagonal, and each lasing mode becomes distributed over many CF states. Sufficiently far 
above threshold, due to these interactions, the lasing modes should lose all resemblance to any 
single CF state. The modes that turn on at higher pump powers tend to be more distributed over 
several CF states even at threshold due to the less diagonal character of their higher threshold 
matrices. This is illustrated in Fig. [S3l which shows the calculated decomposition in CF states at 
Dq/Dqc = 123.5035 of the green mode in Fig. 2. For the mode chosen there are four CF states 
with a weight greater than 10%. The 3rd CF state here has the maximal contribution of 25.6%, 
while it is 55.8% at Dq/Dqc = Dl,^ = 71.9519. The detailed behavior here merits further study 
as to its dependence on size and type of disorder and on kaR. 

A second important observation about the spatial structure of the lasing modes is that the 
amplitude of the lasing modes increase quasi-exponentially towards the gain boundary, see 
Fig. 4. This is not due to high-Q passive cavity modes spatially localized at the boundary (as 
it is for whispering gallery modes of uniform spheres or cylinders); here there are no high-Q 
passive cavity modes at all (see Fig. 1). This effect is due to the high gain needed to initiate 
lasing in such leaky systems. Full understanding of what determines the growth rate in this 
two-dimensional case will also require further study. 

Materials and Methods 5: "Frequency repulsion" in DRLs 

We can now be more precise about the spatial correlation of modes with nearly degenerate 
frequencies. By construction, 7^„(A;) has a real eigenvalue equal to 1/Dq at the frequencies 
kfj_ for each mode that is lasing. If there were two modes with k^ = k^, then the complex random 
matrix Tmn{k) would have an exact accidental degeneracy. Such matrices are a set of measure 
zero in the ensemble, and instead, as the eigenvalues A^, \y approach each other in the complex 
plane, there is an avoided crossing and strong mixing of the eigenvectors a'', a*^ (we see this 
mixing numerically). Strong overlap of a'^, in turn implies strong spatial correlation, strong 
hole-burning interaction and suppression of the weaker intensity mode by the higher intensity 
mode. The net effect is an apparent repulsion between lasing frequencies, even though they 
are not themselves the eigenvalues of a random matrix. The crossing of lasing and non-lasing 
frequencies (see Fig. 2) does not require degeneracy of the threshold matrix and is allowed. 
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Figure S 1 : Eigenvalue flow for the non-interacting threshold matrix. Colored circles represent 
the trajectories in the complex plane of \/\^{k) as a function of k for a scan k = {ka — 
7±, ka + 7±) {kaR = 30,7_Li? = 1). Shown are only a few eigenvalues out of Ncf = 16 
modes included. The direction of flow as k is increased is indicated by the arrow. Star-shapes 
mark the initial values of 1/X^{k) at = ka — 'y±. The full red circles mark the values at 
which the trajectories intersect the real axis, each at a different value of k = k^, providing the 
non-interacting thresholds d\'^ . 
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Figure S2: Evolution of the thresholds and lasing frequencies as a function of Dq. (A) Evolution 
of the quantities X^f!^''Do{i) as the pump is increased. A^t''' denotes the real calculated at step 
i of the discretization of the full Do-range. At a given step i, all modes ^ which are below the 
threshold (delineated by the line A^Dq = 1) are non-lasing at pump Do(i). Once a mode starts 
to lase its corresponding value X^Dq is clamped at X^Dq = 1. The full colored lines represent 
the modes which start lasing within the calculated range of Dq. The dashed lines represent the 
evolution X^/^^DQ{i) i.e. the evolution, had mode-mode interactions not existed, while gray lines 
indicate modes which do not turn on in the pump range shown. The color-coding is identical 
to that of Fig. 2. (B) Analogous evolution of lasing frequencies at which X^^\k) become 
real. The non-lasing modes are drawn in dashed lines, turning into full lines as the modes begin 
to lase. Gray lines represent modes which never lase in the calculated range of Dq. Note the 
absence of frequency repulsion of non-lasing modes. 
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CF3;2S.6% CF6;19.0% CF 11; 10,4% CF 14; 6.5% 



Figure S3: Spatial structure of lasing modes. False color plot of the spatial structure of 
the third mode (green curve in Fig. 2) and its expansion in the CF basis (upper right) at 
Dq/Dqc = 123.5035. The structures of four of the CF states with largest weights are shown 
at the bottom. The maximum amplitudes of the field distributions in these plots are rescaled to 
unity to facilitate comparison. 
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